The ReLPM Exponential Integrator for FE Discretizations of Advection-Diffusion Equations

We implement an exponential integrator for large and sparse systems of ODEs, generated by FE (Finite Element) discretization with mass-lumping of advection-diffusion equations. The relevant exponentiallike matrix function is approximated by polynomial interpolation, at a sequence of real Leja points related to the spectrum of the FE matrix (ReLPM, Real Leja Points Method). Application to 2D and 3D advection-dispersion models shows speed-ups of one order of magnitude with respect to a classical variable step-size Crank-Nicolson solver. 1 The Advection-Diffusion Model We consider the classical evolutionary advection-diffusion problem    ∂c ∂t = div(D∇c)− div(cv) + φ x ∈ Ω, t > 0 c(x, 0) = c0(x) x ∈ Ω c(x, t) = gD(x, t) x ∈ ΓD, t > 0 〈D∇c(x, t), ν〉 = gN(x, t) x ∈ ΓN, t > 0 (1) with mixed Dirichlet and Neumann boundary conditions on ΓD ∪ ΓN = ∂Ω, Ω ⊂ IR, d = 2, 3; cf. [1]. Equation (1) represents, e.g., a simplified model for solute transport in groundwater flow (advection-dispersion), where c is the solute concentration, D the hydrodynamic dispersion tensor, Dij = αT|v|δij + (αL − αT)vivj/|v|, 1 ≤ i, j ≤ d, v the average linear velocity of groundwater flow and ⋆ Work supported by the research project CPDA028291 “Efficient approximation methods for nonlocal discrete transforms” of the University of Padova, and by the subproject “Approximation of matrix functions in the numerical solution of differential equations” (co-ordinator M. Vianello, University of Padova) of the MIUR PRIN 2003 project “Dynamical systems on matrix manifolds: numerical methods and applications” (co-ordinator L. Lopez, University of Bari). Thanks also to the numerical analysis group at the Dept. of Math. Methods and Models for Appl. Sciences of the University of Padova, for having provided FE matrices for our numerical tests, and for the use of their computing resources. M. Bubak et al. (Eds.): ICCS 2004, LNCS 3039, pp. 434–442, 2004. c © Springer-Verlag Berlin Heidelberg 2004 The ReLPM Exponential Integrator 435 φ the source (cf. [2,3]). The standard Galerkin FE discretization of (1) nodes {xi} N i=1 and linear basis functions {φi} N i=1 [4], gives the linear system of ODEs { P ż = Hz + b+ q, t > 0